################################################################################
#       Copyright (C) 2010 Michael Yurko <myurko@gmail.com>
#
#This file is part of pyQAP.
#
#pyQAP is free software: you can redistribute it and/or modify
#it under the terms of the GNU General Public License as published by
#the Free Software Foundation, either version 2 of the License, or
#(at your option) any later version.
#
#pyQAP is distributed in the hope that it will be useful,
#but WITHOUT ANY WARRANTY; without even the implied warranty of
#MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#GNU General Public License for more details.
#
#You should have received a copy of the GNU General Public License
#along with pyQAP.  If not, see <http://www.gnu.org/licenses/>.
################################################################################

cimport numpy as np
cimport cython

ctypedef Py_ssize_t index_t
ctypedef Py_ssize_t np_int_t
ctypedef Py_ssize_t int_t
ctypedef np.float64_t float_t
ctypedef Py_ssize_t bool_t

import numpy as np
import classes
import lap
import cython
'''
Contains helper functions for bounds.
'''

@cython.boundscheck(False)
@cython.wraparound(False)
def gilmore_lawler_full(problem):
    '''
    Return the Gilmore-Lawler lower bound for the given f and d matrices. Uses
    LAPJV to solve the produced LAP.
    
    Examples:
    
    >>> from pyQAP import *
    >>> lap.create_hungarian(100)
    >>> p = problems.QAPLIB('had16')
    >>> from bound_helpers import gilmore_lawler_full
    >>> gilmore_lawler_full(p)
    3358.0
    >>> p = problems.QAPLIB('kra30a')
    >>> gilmore_lawler_full(p)
    68360.0
    >>> p = problems.QAPLIB('nug30')
    >>> gilmore_lawler_full(p)
    4539.0
    >>> p = problems.QAPLIB('rou20')
    >>> gilmore_lawler_full(p)
    599948.0
    >>> p = problems.QAPLIB('tho30')
    >>> gilmore_lawler_full(p)
    90578.0
    
    '''
    cdef np.ndarray[float_t, ndim=2, mode='c'] f = problem.f
    cdef np.ndarray[float_t, ndim=2, mode='c'] d = problem.d
    cdef np.ndarray[float_t, ndim=2, mode='c'] c = problem.c
    cdef index_t n = problem.n
    cdef index_t i, j, k, l, m
    cdef float_t s
    #get sorted version of flow and distance matrices
    cdef np.ndarray[float_t, ndim=2, mode='c'] fn = np.empty((n,n-1))
    cdef np.ndarray[float_t, ndim=2, mode='c'] dn = np.empty((n,n-1))
    for i in range(n):
        k = 0
        for j in range(n):
            if i != j:
                fn[i,k]=f[i,j]
                dn[i,k]=d[i,j]
                k+=1

    fn.sort(axis=1)
    dn.sort(axis=1)
    
    #create cost matrix
    cdef np.ndarray[float_t, ndim=2, mode='c'] cost_matrix = np.empty((n, n))
    for i in range(n):
        for j in range(n):
            #add initial
            s = f[i, i]*d[j, j]
            #setup for indices
            k = 0
            l = n-2
            f_i = range(n-1)
            d_j = range(n-2, -1, -1)
            #add innerproduct
            for m in range(n-1):
                s += fn[i, k]*dn[j, l]
                k += 1
                l -= 1
            #add linear cost
            s+=c[i, j]
            #assign cost
            cost_matrix[i, j]=s

    #compute LAP cost and return
    return lap.hungarian(cost_matrix,  n)

def delete_ij(m, i, j):
    '''
    Returns a view of the matrix without row i and colum]n j. Note: this does not
    make a copy. If the resulting array is changed, then the original array will
    also be changed.
    '''
    m_view = m[:-1,:-1]
    m_view[:i,j:] = m[:i,j+1:]
    m_view[i:,:j] = m[i+1:,:j]
    m_view[i:,j:] = m[i+1:,j+1:]
    return m_view

@cython.boundscheck(False)
@cython.wraparound(False)
def sub_problem(original, psol):
    '''
    Generates a sub-problem from the original given a partial solution.
    '''
    cdef index_t i, j, z, q, r, x, placed, n
    cdef float_t offset
    placed = psol.placed
    cdef np.ndarray[int_t, ndim=1, mode='c'] sol = psol.sol
    cdef np.ndarray[float_t, ndim=2] c_scratch =\
                                              np.empty((original.n, original.n))
    cdef np.ndarray[float_t, ndim=2] f
    cdef np.ndarray[float_t, ndim=2] d
    cdef np.ndarray[float_t, ndim=2] c
    
    f = original.f.copy()
    d = original.d.copy()
    c = original.c.copy()
    n = original.n
    
    #the intermediate reduced problem
    #keep track of actual indices
    cdef np.ndarray[int_t, ndim=1, mode='c'] reduced = np.zeros((original.n),
                                                                 dtype=np.int)
    for z in range(placed):
        #get row and column to delete
        q = 0           #since working off of already reduced matrix
        r = -1
        for x in xrange(sol[z]+1):
            if reduced[x]==0:
                r+=1
        reduced[sol[z]] = 1        #mark as reduced
        #compute new c matrix
        offset = c[q,r]/(n-1)
        for i in range(n):
            for j in range(n):
                c_scratch[i, j] = offset \
                +c[i, j]\
                +f[i, q]*d[j, r]\
                +f[q, i]*d[r, j]
        #remove given row and column
        f = delete_ij(f, q, q)
        d = delete_ij(d, r, r)
        c_reduced = delete_ij(c_scratch, q, r)
        c = c_reduced
        n -= 1
        
    return classes.Problem(f.copy(), d.copy(), c.copy())
        
def anstreicher_brixius(f, d, c, n):
    #create basis
    v = np.asmatrix(np.empty((n, n-1)))
    v[1:, :] = np.identity(n-1)
    v[0, :] = -np.ones((1, n-1))
    
    #create a_hat and b_hat
    a_hat = v.T*f*v
    b_hat = v.T*d*v
    
    #get eigenvale decompositions
    lambda_ = np.linalg.eigvals(a_hat).real
    sigma = np.linalg.eigvals(b_hat).real
    
    #sort in opposite orders
    sigma.sort()
    lambda_.sort()
    lambda_ = lambda_[::-1]     #reverse array
    
    #create LAP matrix
    c = np.empty((n-1,n-1))
    for i in xrange(n-1):
        for j in xrange(n-1):
            c[i,j] = lambda_[i]*sigma[j]
            
    #create constraints for LAP as LP
    a_lap = np.zeros((2*n, n*n))
    for i in xrange(n):
        for j in xrange(n):
            a_lap[i, (i-1)*n+j] = 1
            a_lap[n+i, (j-1)*n+i] = 1
    
    b_lap = np.ones((2*n,1)) 
